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Universal energy diffusion in a quivering billiard 
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We introduce and study a model of time-dependent billiard systems with billiard boundaries 
undergoing infinitesimal wiggling motions. The so-called quivering billiard is simple to simulate, 
straightforward to analyze, and is a faithful representation of time-dependent billiards in the limit 
of small boundary displacements. We assert that when a billiard’s wall motion approaches the quiv¬ 
ering motion, deterministic particle dynamics become inherently stochastic. Particle ensembles in 
a quivering billiard are shown to evolve to a universal energy distribution through an energy diffu¬ 
sion process, regardless of the billiard’s shape or dimensionality, and as a consequence universally 
display Fermi acceleration. Our model resolves a known discrepancy between the one-dimensional 
Fermi-Ulam model and the simplified static wall approximation. We argue that the quivering limit 
is the true fixed wall limit of the Fermi-Ulam model. 

PACS numbers: 05.45.-a, 05.40.-a 
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I. INTRODUCTION 

Billiards are remarkably useful physical models; they allow a diverse range of classical dynamics to be understood 
intuitively through easy-to-visualize particle trajectories and are a natural setting for quantum and wave chaos [I|, 
while the discrete time nature of particle-billiard boundary interactions make classical billiards especially amenable to 
numerical study. Time-dependent billiards (billiards with boundaries in motion) in particular can be found in a wide 
range of applications: KAM theory f^-H . one-body dissipation in nuclear dynamics [g, Fermi acceleration [3. l6l-[T^ . 

and adiabatic energy diffusion [3, , for example. 

This work was originally motivated by the desire to study and simulate classical particle trajectories in time- 
dependent billiard systems. The task is complicated by the boundary’s displacement, which produces implicit equa¬ 
tions for the time between particle-boundary collisions. We propose a fixed wall simplification by considering the 
limit of infinitesimally small boundary displacements. Our limit will be called the quivering limit, and the resulting 
billiard system will be called a quivering billiard. The purpose of this paper is to show that, although simple, quiver¬ 
ing billiards are accurate descriptions of time-dependent billiards in the limit of small boundary displacements, and 
to support our conjecture that any physically consistent, non-trivial, fixed wall simplification of a time-dependent 
billiard must be physically equivalent to a quivering billiard. Using physical reasoning, we will argue that in the 
quivering limit, deterministic billiard dynamics become inherently stochastic. Then, by utilizing the simplifications 
allowed by stochastic methods and fixed billiard walls, we will derive analytic expressions to describe energy evolution 
in a quivering billiard. Our investigations will uncover universal behavior in time-dependent billiards when billiard 
motion is close to the quivering limit, and our results will enable us to addresses several issues that have been raised 
in previous Fermi acceleration and time-dependent billiard literature. 

The outline of this paper is as follows. In Sec. m we first define a quivering billiard and determine its behavior in one 
dimension, and then generalize to quivering billiards in arbitrary dimensions. The energy statistics of a single particle 
and a particle ensemble are examined in Sec. Ell and the results are discussed in the context previous literature in 
Sec. lIVI In Sec.|Vl we give examples of quivering billiards and present numerical analyses, and we conclude in Sec. lVIl 


II. THE QUIVERING LIMIT 

In this section, we define quivering as a particular limit of time-dependent billiard motion. Because the dynamics 
are so poorly behaved in this limit, billiard systems can only be described stochastically. For simplicity, we first work 
with a one-dimensional billiard with a single moving wall, and then extend to arbitrary billiard motion in arbitrary 
dimensions. 


A. The 1-D Fermi-Ulam Model 


Consider a particle in one dimension bouncing between two infinitely massive walls. One wall is fixed at a; = 0, 
and the other oscillates about its mean position ai x = L, where we take L > 0. The particle’s energy fluctuates 
due to collisions with the moving walk and the dynamical system corresponding to the particle’s motion defines 
the well-known Fermi-Ulam model Suppose that the moving wall oscillates periodically with period r, 

characteristic oscillation amplitude a, and characteristic speed Uc = a/r. The moving wall’s position x(t) and velocity 
u(t) at time t can be written as 


x{t) 

u{t) 


L -I- g{t) 
dt ’ 


( 1 ) 


where g{t) is some piecewise smooth r-periodic function with mean zero. The wall velocity scales like Uc, and g{t) 
scales like a. To make the scaling obvious, we note that g{t) depends on t only through the value of t mod r, and we 
make the following substitutions: 


T(t) = - mod 1 (2) 

T 

g{t) = ah{^{t)). 


The quantity 4'(t) will be referred to as the wall’s phase. Here, h is regarded as a function of di', and means 

evaluated for = 'l'(t). The quantity is just g{t) rescaled to have a characteristic oscillation amplitude 
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FIG. 1. Spacetime diagram over one period of a moving wall’s motion. The smooth curve represents the wall’s position, and 
particles approach the wall along the diagonal arrows to collide at times ta,ta + and ti, + St 


of unity. The state of the wall at time t is thus 


x{t) = L + a (3) 

u{t) = Uc h'(T'(t)), 

where the h' denotes the derivative of h with respect to its argument 'll. 

We define the quivering limit of the Fermi-Ulam model by taking a, r —>■ 0 while holding Uc constant and leaving 
the dependence of h on th fixed. In the quivering limit, the moving wall’s position reduces to x{t) = L, so no implicit 
equations for the time between collisions arise from the dynamics. This simplification comes at a price; when r —>■ 0, 
d' oscillates infinitely fast in time, and u{t) does not converge to any value for any given t. That is, in the quivering 
limit, u{t) becomes ambiguous to evaluate. Our task now is to physically interpret and resolve this ambiguity. 

Note that in the quivering limit, the wall makes infinitely erratic motions at finite speeds; the derivative of 
g{t), scaling like ajr'^, diverges for all n > 2. An infinitesimal change in the state of a particle results in a finite 
and essentially unpredictable change in the wall’s velocity at the time of the next bounce. We assert that one could 
never, even in principle, specify the state of the particle with enough precision to reliably predict the velocity of the 
moving wall, and thus the change in particle energy, during the next collision. We therefore claim that in the quivering 
limit, the dynamics of the Fermi-Ulam model become inherently stochastic; deterministic particle trajectories defined 
on phase space transition to stochastic processes defined on a probability space. Given any initial condition, the 
resulting particle trajectory actually represents one possible realization drawn from an ensemble of initial conditions 
inhnitesimally displaced from one another. The wall’s velocity during a collision will be treated as a random variable, 
and we now hnd the corresponding probability distribution. 

Consider again the moving wall with non-zero a and t. Let P(u\0) be the probability density for a stationary 
observer to measure the velocity u during a randomly timed snapshot of the wall: 

PfulO) = — [ dtS(u — u(t)) (4) 

T Jo 

1 

d'i>6{u — Uc /i'(4>)). 

The reason for placing the conditional |0 in the argument of P will become apparent shortly. We note that P(u|0) is 
normalized, so it is indeed a well-defined probability density. In the quivering limit, Uc and the dependence of h on 41 
remain constant, so P{u\0) remains well-defined and unchanged. If the stationary observer were to measure the wall 
velocity in the quivering limit, any observation, no matter how well-timed, would be an essentially random snapshot 
due to the wall’s infinitely erratic motion. We thus take P(u|0) to be the probability for a stationary observer to 
measure the wall with velocity u when the wall is quivering. 

The particle bouncing between the walls effectively measures the wall’s velocity during collisions, but the particle 
is not a stationary observer. Collisions with large relative speeds of approach occur more frequently than collisions 
with small relative speeds of approach, so there exists a statistical bias that favors collisions for which the wall moves 
towards the particle. If the quivering dynamics are to be physically consistent with the Fermi-Ulam dynamics, this 
statistical bias must be incorporated into the probability distribution used to determine the wall’s velocity during 
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collisions. The mathematical realization of the statistical bias can be found with the aid of Fig. [U a construction first 
employed by Hammersley Q and Brahic 0. 

In Fig. [ll the position of the moving wall in the Fermi-Ulam model is plotted over one period of motion in the 
interval {to, to + t). Consider an ensemble of particles approaching the moving wall with speed v. For the moment, 
we assume that v is larger than the maximum wall velocity Umax- The particles are launched from x = 0 at a uniform 
rate over a period of duration r such that they all collide with the wall during the interval {to, to -|- r). We concern 
ourselves only with the first collision each particle makes with the moving wall. Four trajectories from the ensemble 
are shown in Fig. [TJ representing collisions with the wall at times ta, ta + St, tb, and tb + St. Because the launch 
times are uniformly distributed, the fraction of particles that collide with the wall between ta and ta St will be 
proportional to the interval Sta = St — Aa- Likewise, the fraction that collide between tb and tb + St will be proportional 
to Stb = St + Ab- Using the geometry of Fig. [T] and the fact that tan(6*) = v, we hnd the probability density for 
randomly selected ensemble member collide with the moving wall at a time t within the interval {to, to + t) to be 

P(u(t)|u) = i(^l-^). (5) 


Multiplying by a delta function and integrating Eq. ([5]) over a period of the wall’s motion gives P{u\v), the probability 
density for a randomly selected ensemble member’s collision to occur when the wall moves with velocity u: 


P{u\v) 


T 

— J dt S{u — u{t)) ^1 
0 

1 

J d'd> S{u — Uc h'{d/)) 
0 

P{u\0) 



Uch'{^) S^ 


( 6 ) 


Because the wall’s average displacement over one period of motion is zero, the product uP{u\0) integrated over all 
wall velocities must also give zero, and P{u\v) is therefore normalized and a well-defined probability density. The 
distribution P{u\v) has a statistical bias towards larger negative u due to the flux factor 1 — u/v. We will henceforth 
refer to P{u\0) as the unbiased distribution and P{u\v) as the biased distribution. In the quivering limit, P{u\v) 
remains well-defined and unchanged. As r —^ 0, an ensemble of particles launched over a period of wall motion from a 
fixed X is essentially equivalent to an ensemble of infinitesimally displaced initial conditions. We therefore take P{u\v) 
to be the conditional probability density to observe a quivering wall with velocity u during a collision, given that the 
particle approaches the wall with speed v > Umax- 

If a particle approaches the moving wall with speed v < Umax-, then P{u\v) will become negative for some values 
of u, and Eq. dH]) will make no sense as a probability density. These u values correspond to impossible collisions for 
which the wall moves with positive velocity away from the particle faster than the particle moves toward the wall. 
Such collisions occur with probability zero, and we can account for this by simply attaching a step-function to the 
biased distribution, yielding 


P{u\v) 


P(u|0)(l-^), v>u, 

N{v)P{u\0) (l — 0(u — u), V < u, 


( 7 ) 


where 0(x) is the unit step function (equal to 0 for x < 0 and 1 for x > 0) and N{v) is a u dependent normalization. 

Equation © determines the statistics of a particle’s energy evolution in a quivering Fermi-Ulam system. As with 
any billiard system, the particle’s energy is simply the kinetic energy where m is the particle’s mass and v is 

its speed. The particle bounces between the two walls as if the system were time-independent, but when colliding 
with the quivering wall at an incoming speed Vi (the particle moves in the positive x direction to collide with the 
moving wall, so Vi is also the incoming velocity), a value for the wall velocity u is selected using the biased distribution 
P{u\vi). The particle’s velocity just after the collision, Vf, is given by 

Vf=2u-Vi, ( 8 ) 

and the corresponding energy change, AE, is given by 


AE = 2mu^ — 2muvi. 


( 9 ) 
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Equations ([8|) and ([9]) are determined using the standard collision kinematics for a particle in one-dimension colliding 
elastically with an infinitely massive moving object. 

Before moving on to higher dimensions, we must address the possibility of particles escaping the billiard interior. 
This issue will plague any fixed wall simplification of time-dependent billiards, and is discussed in detail in Ref. [l^ . 
From Eq. ®, we see that if 0 < it < < 2u, the particle does not turn around after a collision with the moving 

wall, but instead slows down and continues forward. We refer to these types of collisions as glancing collisions. For 
non-zero a and r, just after a glancing collision, the particle continues forward slower than the wall moves outward, 
so the particle will remain within the billiard interior. With a fixed wall simplification, however, the wall does not 
actually move outward after a glancing collision, so the particle will continue forward and escape the billiard interior. 
A particle escaping through a hard wall is a non-physical by-product of setting a = 0, so in order to make a physically 
reasonable fixed wall simplification, one must always devise a method to handle glancing collisions. Our method for 
a quivering Fermi-Ulam system is devised as follows. 

For non-zero a and t, after a glancing collision occurs, the wall continues to evolve through its period, and one of 
two possibilities will occur. The wall may slow down sometime after the glancing collision and allow the particle to 
catch up and collide again, or the wall may reverse its direction and move inward sometime after the glancing collision, 
also allowing the particle to collide again. In either case, a second collision occurs after the first collision, and as a and 
T approach zero, the second occurs essentially instantaneously after the first. Therefore, we treat a glancing collision 
in a quivering Fermi-Ulam billiard as a double collision. When a particle with speed Vi (also the particle’s velocity) 
collides with the quivering wall, we draw a u value from the distribution P{u\vi). If the selected value of u is such 
that 0 < u < Vi < 2u, the particle’s new speed Vf (also velocity) is given hy Vf = 2u — Vi, and we draw a new u value 
from the distribution P{u\vf). If the second u value gives another glancing collision, we again update the particle’s 
speed and then draw a third u value. The process is repeated until a non-glancing collision occurs, and the whole 
event (which occurs instantaneously) is treated as a single collision. 


B. Arbitrary Time-Dependent Billiards 

We now generalize to arbitrary billiards in arbitrary dimensions. Consider a time-dependent billiard in d dimensions 
moving periodically through some continuous sequence of shapes with period r, characteristic oscillation amplitude a, 
and characteristic speed Uc = a/r. The evolution of any one point on the boundary will be denoted by the path q(t), 
where q(t -I- r) = q(<). For every t, the set of all boundary points {q(t)} is assumed to define a collection of unbroken 
d — \ dimensional surfaces, which we refer to as the boundary components, enclosing some d dimensional bounded 
connected volume. The outward unit normal to the billiard boundary at the point q(t) is denoted by n(q(t)), and 
the velocity of the boundary point q(t) is denoted by u(q(t)) = dq(t)/dt. The billiard shape evolves continuously in 
time, and we assume that the boundary components remain unbroken throughout their evolution, so u(q(t)) forms 
a smooth vector field with domain on the boundary {q(t)} for any fixed t. Likewise, n(q(t)) forms a smooth field 
on {q(t)} for any fixed t, except possibly at corners, where n(q(t)) is ill-defined and discontinuous. We denote the 
outward normal velocity of the point q(t) by u(q(t)) = u(q(t)) • n(q(t)). 

Denote by q the average of q(t) over one period: 


^=- f dtq(t). (10) 

''' Jo 

Noting that the boundary components remain unbroken throughout the period of motion, it is straightforward to show 
that set of average boundary points {q} forms a collection of unbroken d— 1 dimensional surfaces. The trajectory q(t) 
and normal velocity it(q(t)) of any given boundary point can be written as functions of the corresponding average 
location q and the time t: 


q(t) = q-bg(q,t) (II) 

u(q,t) = 9tg(q,i) • n(q(t)), 

where g(q, t) is a piecewise smooth in time r periodic function with a time average of zero. g(q, t) scales like a and 
u(q(t)) scales like Uc- Equation (ITT|) depends on t only through the value of 4'(t) = t/r mod I, so we write 

q(t) = q-bah(q,^'(t)) (12) 

u(q, t) = Uc 9^h(q, 5'(t)) • n(q(<)). 

where ah(q, ^'(t)) = g(q, t). Analogously to the one dimensional case, h is regarded as a function of q and 'k, and 
h(q, 5'(t)) means h(q, 'k) evaluated for = 'l'(t). The quivering limit of an arbitrary dimensional billiard is defined 
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FIG. 2. Collision geometry in a two-dimensional billiard. A particle with velocity v approaches the point q on the billiard 
boundary, where the outward unit normal vector is n. The dotted line represents the tangent line to the boundary at q 


by taking a,T — ^ 0 while holding Uc and the dependence of h on T' and q constant. In this limit, the billiard’s 
boundary points become fixed in time at the average locations {q}, so the outward normal vectors become fixed in 
time as well. Thus, in the quivering limit, we have 

q(t) = q (13) 

■u(q,t) = ■Uc5^h(q, T'(<)) • n(q) 

= Uc /i'(q,T'(t)), 

where we write h'{q, T'(t)) = (9^h(q, T'(t)) • n(q) for brevity. Any time-dependent billiard taken to the quivering limit 
will be called a quivering billiard. 

Analogously to the one dimensional case, we define the unbiased distribution for each q: 


P(u|0,q) 


dT'^(u 


Uc /i'(q, T')). 


(14) 


The biased distribution for each q can also be defined analogously to the one dimensional case, but we must also 
consider the collision angle 6, depicted for two-dimensional billiard in Fig. [51 For a particle approaching the boundary 
point q with speed v, 0 is the angle between the particle’s velocity vector and the d — 1 dimensional tangent surface 
to the wall at q, and v sin(0) thus gives the component of the particle’s velocity in the n(q) direction. If the particle 
collides when the wall has normal velocity m, then the relative speed of approach just before the collision is determined 
by V sin(0) and m, so v sin(0) determines the statistical bias towards collisions with large negative u. We account for 
this by simply replacing v with z;sin(0) in Eq. ([7]), yielding 


P{u\v,(i, 0) 


P{u\0,q) (l - , vsin{0) > Umaxin) 

N{v,0)P{u\O,q) (l - 0(wsin(6») - u), vsin{0) < Umax{(l), 


(15) 


Equation m determines the statistics of a particle’s energy evolution in a quivering billiard. 

To summarize, we describe how one may construct a quivering billiard and determine a particle’s trajectory, without 
the need to define a real, fully time-dependent billiard and take the quivering limit. First, one must select a billiard 
shape by defining a surface {q}, then set boundary quivering by giving a value to Uc and defining a scalar field 
/i'(q, 4') on {q}. If the constructed quivering billiard is to honestly represent some deterministic billiard’s motion in 
the quivering limit, then h'{q, d^) should be chosen to be a smooth function of q for any ik wherever n(q) in continuous. 
Using the field h' and the value of Uc, one may then calculate the unbiased distribution P(u|0,q) from Eq. (fHl) for 
any q on the billiard boundary. For a particle in free flight inside the quivering billiard, the next collision location is 
found deterministically using the geometry of the billiard boundary, just as with a time-independent billiard. When 
a particle with velocity and speed Vi collides with the boundary at q with a collision angle 0i, we draw a value of 
u from the distribution P(M|ui, q, 0i). The particle’s velocity component tangent to the boundary remains constant, 
and the component normal to the boundary just after the collision, v/ • n(q), is given by 


V/ • n(q) = 2u-Vi- h(q) 
= 2u — Vi sin(0i). 


(16) 




7 


The corresponding change in energy, Ai?, is given by 

AE = 2mu^ — 2mu Vi sin {6i). (17) 

Analogously to the one dimensional case, if the selected value of u is such that 0 < u < Vi sin(0i) < 2 m, then a glancing 
collision occurs, and we draw a second value of u using the same collision located and updated particle speed and 
collision angle, determined from Eqs. (HH) and dnl)- 


III. ENERGY STATISTICS 

In this section, we study in detail the statistical behavior of particles and ensembles in a d-dimensional quivering 
billiard, with the aim of describing energy evolution of a ensemble of initial conditions as a diffusion process. Our 
notation will be as follows: q?, is the location of a particle’s collision with the billiard boundary, 9b is the 
collision angle, Ub is the selected value of the wall velocity during the 6*^ collision (sampled using Eq. (ITKll b vt-i is 
the particle’s speed just before the collision, and AEb is the change in particle energy due to the 6*^ collision, 
given by 


AEb = ^mul — 2mub Vb-i sin {9b) ■ (18) 

In order to derive analytic results, we will assume that the initial particle speeds vg are much larger than Uc, and 
we will solve to leading order in the small parameter e = Uc/vg. We regard Uc as an 0(1) quantity, and vg as an 
0{£~^) quantity. This approximation allows us to ignore glancing collisions in our analysis, and also allows us ignore 
the possibility of Vb-i sin{9b) < Umax{(lb), so that the biased distributions at the time of collision always take the form 
P{ub\vb-i,qb, 9b) = P{ub\0, Qb) (1 — Ubjvb-i sin(db)) (as opposed to the more complicated Eq. (fT^ b The assumption 
e <C 1 is not particularly restrictive; even if particles begin with an initial speed comparable to or less than Uc, energy 
gaining collisions are more likely than energy losing collisions due to the flux factor in the biased distribution, and 
a slow particle will gain roughly mu^ of energy during a collision according to Eq. m- Therefore, a slow particle 
will more than likely gain speed Uc 0(1) during a single bounce, and after 1/S bounces, where <5 <C 1 is some small 
number, the particle will more than likely have a speed v such that Uc/v <C 1. Thus, slow particles are very likely 
to eventually become fast particles, and the assumption Uc/v <C 1 will give a better and better approximation over 
time. 

In the analysis, it will prove useful to consider both the full dynamics and frozen dynamics, as is done in Refs. [l^ . [l^ . 
If the frozen dynamics are used at the collision, the energy change AEb is calculated, but the particle’s energy 
remains constant, and the angle of reflection is equal to the collision angle 9b- In other words, the frozen dynamics 
are identical to those of a time-independent billiard, but we calculate and keep track of the AEbS that would have 
occurred had the billiard walls been quivering. In the full dynamics, the particle’s energy is actually incremented by 
the calculated value of AEb, and the angle of reflection is consequently altered. 


A. Expectations 


Consider single a particle with energy Eg released at time tg in a d-dimensional quivering billiard. The resulting 
particle trajectory generates a sequence of energy increments {AEi, AE 2 ,..., AEb-i, AEb, AEb+i, ...}. Let the opera¬ 
tor {...}b denote the conditional expectation value of the quantity ..., given the outcomes of the previous 6—1 bounces. 
The first 5—1 bounces determine Vb-i, qb, and 9b, so the 5‘^ conditional expected energy change, fib = {AEb}b, can 
be calculated using the biased distribution P{ub\vb-i,qb,db) and the expression for AEb in Eq. (fTSl) : 


9b — {Aii/b}b 


(19) 


= J dubP{ub\vb-i,qb,db)AEb 

= / dubP{ub\0,qb) i dmul - ■ . 

J \ Vb-isin{ 9 b) 


2mul 


— 2mub Vb-i sin(db) 


The integral in Eq. (fT^ is taken over all possible values of Ub at qb. 

Let M„(qb) denote the moment of the wall velocity at qb as measured by a stationary observer: 


M„(q) = 


dMP(M|0, q)M". 


( 20 ) 



By construction, Mi(q) = 0 for all q. Otherwise, Mn{c[b) generally scales like u'^. The conditional mean thus simplifies 
to 


tJ-b= dubP{ub\0,q_b)4:mul ( 1 - --^ 

J \ 2vb-i sin(6>b) 


= 4mM2(qf,) 1 - 


M3(qb)/M2(qb) 

2vb-i sm{9b) 

Similarly, the conditional variance is given by 

^ {{AEbr}b - {AEb}l 

= J dubP{ub\vb-i,qb,db) {{AEbf - {AEb}l) 

Vb-i sin{9b) 


m 2 I sin (fi'fc) 


= 4TO^[M2(q6)] 


- 3t 


+ 3' 


Mi{qb) 


-4- 


y M2(qb) [M2(qb)]2/M3(qb) [M2(qb)]2 

4M3(qb)/M2(qb) - M5(qb)/[M2(qb)]2 [M3(qb)]7[^2(qb)]'' 


Wb-i sin(6»b) 


7-1 sin (^b) 


( 21 ) 


( 22 ) 


The terms enclosed in the parentheses of Eqs. (f?!]) and (l2^ are ordered in increasing powers of e. To leading order, 
we have 


fib = 4mM2(qb) (23) 

al = 4TO^M2(qb)7-iSi7(6'b) 

The quantities ^b and are 0(1) and O (e~^), respectively; average energy gain is moderate, and fluctuations are 
huge. 


B. Correlations 


The conditional covariance between adjacent bounces, Covb,b+i, is defined by 

Covb,b+i = {{AEb - {AEbIb) {AEb+i - {A£;b+i}b)}b (24) 

= {AEbAEb+i}b — {AEb}b{AEb+i}b- 

The conditional expectations in Eq. (1241) are taken given the outcomes of the previous b— 1 collisions, with the outcome 
of the collision yet to be determined. That is, we must average over all possible realizations of the stochastic process 

Eb-i —>■ Eb-i+ AEb Eb-i+AEb + AEb+i, given the first b— 1 collisions. Denote {AEb+i\ub}b+i as the conditional 
expectation of Eb+i, given the first 6—1 collision outcomes and supposing that Ub is the wall velocity during the 
collision. The expression for {AEb+i}b is then 

{AEb+i}b = J dubP{ub\vb-i,qb,Ob){AEb+i\ub}b+i- (25) 

The expression for {AEbAEb+i}b can be written similarly: 


{AEbAEb+i}b = J dubdub+iP{ub\vb-i,qb,db)P{ub+i\vb,qb+i,Ob+i\ub)AEbAEb+i (26) 
= J dubP{ub\vb-i,qb, Ob)AEb{AEb+i\ub}b+i- 

The term P{ub+i\vb, qb+i, fi^b+il^b) denotes the value of P{ub+i\vb, qb+i, 0b+i) when Vb, 6b+i, and qb+i are determined 
given the first 6—1 collision outcomes while supposing that Ub is the wall velocity upon the P^ collision. Equation (IMl) 
can thus be expressed as 

Covb,b+i = J dubP{ub\vb-i,qb, 0f)){A£’b+i|Mb}b+i {AEb — {AEb}b) ■ (27) 

If the frozen dynamics are used at the P^ collision, then Vb, Ob+i, and qb+i are independent of Ub, so we have 


{AEb+i\ub}b+l\F = {Ai?b+i}b+l|F = 9b+l\F, 


(28) 
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where ...|f denotes the quantity ... evaluated using the frozen dynamics. ^h+i\F carries no Ub dependence, so it can 
be brought outside of the integral in Eq. (I?71) . giving 


C0Vh,f,+ l|F = 0. 


(29) 


Adjacent energy increments are thus statistically uncorrelated in the frozen dynamics. 

Under the assumption £ <C 1, the frozen dynamics closely resemble the full dynamics over the time scale of a few 
bounces H- Over such a time scale, we can regard the full dynamics trajectory as a stochastic perturbation of the 
deterministic frozen dynamics trajectory. Let qb+ijMf, = q^+ilp + (5q;,_|_i|ub be the {b + 1)*^ collision location when 
the full dynamics are used at the bounce, given the first 6—1 collisions and supposing that Ub is the observed wall 
velocity upon the 6‘^ collision. Equation (l23l) then gives, to leading order in e 

{AEb+i\ub}b+i = ^mM2{qb+i\ub) (30) 

= AmM2{qb+i\F) + 4to VM2(qb+i|F) ■ <5q6+i|u&. 

where the gradient VM 2 is constrained to act along directions tangent to the billiard boundary at qb+ij^. In the 
appendix, we solve for ||(5qt,+i|Mf,|| to leading order in £ and find 


ll<5qf,+i|ub|l 


2L&|f 


cos{9b) |ub| 
sin(0f,+i|F) Vb-i ’ 


(31) 


where Lf,|F is the distance between the 6‘^ and 6+1*^ collision locations in the frozen dynamics. Combining Eqs. (I27L 
dsni), and gives to leading order in e 


Covb,b+i = j dubP{ub\vb-i,(lb,9b) (4to VM 2 (q 6 +i|F) • Sqb+i\ub) (AE& - {AEb}b) 


6q6+i|Mb 


= / dMbP(Mb|0, q&) ( 4mVM2(q&+i|F) , ,, 

X 2£j,|F f^mul - 2m - - 2mubVb-i sm{9b) 

sm[9b+i\F)Vb-iJ \ Vb-isin{9b) 


-4mM2(qb) + 4mM2(qb) 


Ub 


Vb-i sin(6»b) 


1 C 2 r I cos( 6 >b)sin( 6 lb) I 1 Db in 1 

=-16to Lb f— r-775 - p—VM2(qb+i f) • / dwbP(Mb 0, qb)-;^^ - 

sin(6'b+i f) j 6q6+i 


Ubil 


:Ub\Ub\. 


(32) 


All but the leading order terms are dropped in the last line of Eq. (15^ . With exception to the one-dimensional case, 
Covb,b+i is thus an 0(1) quantity. In a one-dimensional billiard, the frozen and full dynamics always give the same 
collision location, so {ALb+i|Mb}b-i-i = {ALb+i}b+i|F) and consequently, Covb^b-i-i is identically zero. 

The conditional correlation Pb,b+i is defined as the normalized conditional covariance, and is given by 


Pb,b+l 


Covb,b-n 

CTbicTb+ljb 


(33) 


To leading order in £, the conditional expectation {ab+i}b can be taken as the frozen dynamics value in Eq. (1331) . 
Therefore, the conditional correlation Pb,b+i is 0{e^) (with exception to the one-dimensional case, where Pb,b+i = 0). 
This quantity is very small, and correlations between more distant collisions will further diminish due to the mixing 
of particle trajectories induced by the stochastic wall motion. We thus conclude that, in any dimension, correlations 
between energy increments effectively decay over the time scale of a single collision. 


C. Ensemble averages 


Consider now a microcanonical ensemble of independent particles with energy Eq released at time tg- The resulting 
trajectories will generate an ensemble of statistically independent energy increment sequences, and we denote AEi^b 
as the 6*^ recorded energy increment of the particle. Define the ensemble averaged 6*^ energy increment (AEb) as 


(AUb) 


AEi^b 

N^oo * ^ N 


(34) 
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and the ensemble averaged conditional mean (/x^) as 


N 

if^b) / ^ ; 

N —¥oo J\l 

2=1 


(35) 


where fii^b = Equation (1^ shows that the conditional variances = {{AEi^b)'^}b — {AEi^bY 

OO 

finite and bounded from above. Noting this, and the fact that the series ^ k~‘^ converges, we deduce 


fc=i 


N 2 

lim ^ < OO. 

AT-J-oo k‘^ 

fe=i 

By Kolmogorov’s strong law of large numbers El , Eq. {|3fij) assures that, with probability unity, 

(AEb) = (fib) ■ 

Combining Eqs. (EZl),®, and ( 1 ^ gives, to leading order in e. 


N 


{AEb) = V 

N^oc ^ 


4mM2(q*,h) 
N 


(36) 


(37) 


(38) 


1=1 

= 4m(M2(qb)) 


where qi^f, denotes the collision location of the particle. By similar law of large number arguments, we also 
have, to leading order in e, 


where 


{ul) = (M2(qb)) 


N 2 

2=1 


and Ui^b is the wall velocity during the 6*^ collision of the particle. To leading order, we thus have 

{AEb) = 4m (ul) . 


(39) 


(40) 


(41) 


D. Energy diffusion 

We now consider the normalized energy distribution of an ensemble of independent particles, denoted by r]{E,t). 
We have thus far shown that energy of any one ensemble member evolves stochastically, in small increments, with 
correlations in energy changes effectively decaying over a characteristic time scale given by time between collisions. 
A particle’s energy evolution is therefore effectively a Markov process describing a random walk along an energy axis, 
so following Refs. [13,Hi! we assert that r]{E,t) evolves like a diffusion process and obeys a Eokker-Planck equation: 

dtr){E,t) = -Oe [giiE,t)r){E,t)] + [g 2 iE,t)r]{E,t)] . (42) 

The functions gi{E,t) and g 2 {E,t), the drift and diffusion terms, respectively, are to be determined in this section. 
The energy of any one particle in a quivering billiard evolves discretely in time, so the continuous time evolution 
implied by Eq. (1421) will be an accurate description of the ensemble only down to a coarse-grained time scale. The 
time scale must be large enough to ensure that most particles in the ensemble experience at least a few bounces off 
the billiard wall, but small enough to ensure the energy change experienced by most particles is small compared to 
their total energy. Generally speaking, a diffusive description of a stochastic process is only accurate over time scales 
larger than the process’s typical correlation time [l3i[I3|- We have established that energy correlations for any one 
particle effectively decay over the time scale of a single collision, thus, the diffusion approach to energy evolution in 
a quivering billiard is justified on any time scale over which g{E, t) can be described by a continuous evolution. 
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The drift term gi{E',t') is defined as the rate of ensemble averaged energy change for an ensemble of particles all 
with energy E' at time t'. Specifically, 




—^0 


{E{t' + At) - Eit')) 
At 


(43) 


where E{t') = E' for all particles in the ensemble, and {E{t' + At)) is the ensemble averaged particle energy at 
time t' + At. We can not actually take the limit At —^ 0 because gi has no meaning over time scales for which the 
evolution of t] appears discontinuous. Instead, we will let At be the average time for which the ensemble members 
make B bounces after time t', and we will find corresponding ensemble averaged change in energy. We assume that 
B is small enough so that the particle energies change very little relative to E' over the time At, so that At is the 
smallest coarse-grained time scale for which Eq. 621) is valid for an ensemble with common energy E'. We let Eb he 
a particle’s energy B bounces after t', and find from Eq. 63 ) 


{Eb -E') = (Y, ^Et 


(44) 


= ^4m(M^) . 

We denote the coarse grained squared wall speed by v?{t'-,B), defined as the time average of (u^) over the first B 
bounces after t'\ 


We thus have 




{Eb - E') = 4m v?{t']B)B 


(45) 


(46) 


The time scale At corresponding to the B bounces after t' is the ensemble averaged total free flight time over which 
the B bounces occur. If we denote by At;, a particle’s free flight time after t', we have 


At = ^ (Atfc). 


(47) 


We are assuming small wall velocities, so the particles’ speeds change very little relative to their initial speed ^2E' jm 
over the B bounces. Therefore, to leading order in e, we have 


Atfc = 


2E' 




(48) 


where Ih denotes a particles 5*^ free flight distance after t'. We now define the coarse grained free flight distance, 
Z(t', B) by time averaging the ensemble average of lb over the first B bounces after t': 




Substituting Eqs. (l4^ and (l4^ into Eq. (l47ll gives 


At = B-l{t'-,B)^l — , 


and substituting for B in Eq. 63 gives 


At 

lit';B) 


(49) 


(50) 


(51) 


Equation (1511) gives the ensemble averaged change in energy over the time At after t' for an ensemble of particles 
with energy E'. Comparing to Eq. (l43l) . we see that dividing both sides of Eq. 63 by At gives us gi{E', t'). We thus 
have, 


^ ^ ^ l{t) 


(52) 
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where we have switched from primed to unprimed variables, and the dependence on B has been suppressed. 
The diffusion term g 2 {E',t') is defined as 


g2iE',t') 


— ^0 


{{E{t' + At)-E{t')f) 
At 


(53) 


where E{t') = E' for all particles in the ensemble, and {E{t' + At)) is the ensemble averaged particle energy at time 
t' + At. An expression for the diffusion term can be found by employing similar methods used to find the drift term. 
Alternatively, 32 (A, t) can be found by invoking Liouville’s theorem, as in Ref. [l^. Combing Liouville’s theorem and 
the Fokker-Planck equation allows one to deduce a fluctuation-dissipation relation: 

= 2^^^ P(i?) 52 (A, t)], (54) 

where S(A) is the microcanonical partition function of a single particle with energy E in the corresponding frozen 
billiard. In a d dimensional billiard, the microcanonical partition function is given by H 

E{E) = ^Vdnd{2m)^ E^-\ (55) 

where ild is the d-dimensional solid angle, and Vd is the d-dimensional billiard’s volume. Combining Eqs. (I511),(IM1), 
and ((55|) . we find 


92{E,t) 


4 4-\/2m u'^(t) ^3 

-r-=- —E^ 

d -|- 1 I (t) 


(56) 


This method of determining g 2 allows for an additive constant, but this constant must be identically zero; when E = 0, 
the particles are motionless and there can be no drift or diffusion of energies, so we must have gi{0,t) = 32 ( 0 , t) = 0 . 
With our expressions for gi and 32 , we may rewrite the Fokker-Planck equation: 


dtv{E,t) 


2a{t) 

d+1 


ds 


l+d / 2-d , 

E—dE [E—g{E,t)j 


(57) 


where we define a{t) as 


_ 4-\/2m u'^{t) 

^ W) 


(58) 


Equation (1571) can be simplified by defining a rescaled time s: 


s = 



(59) 


which gives 


dsiliE, s) 


d+1 


Oe 


1+d / 2-d , 

E—Oe [E—rj{E, s)j 


(60) 


Equation (l60l) can be solved by separation of variables. We assume a solution of the form (j){s)f{E), and upon 
making the substitutions F{E) = E~r- f{E) and z = E^ one finds a first order homogeneous linear differential 
equation for (j){s) and a Bessel equation of order d — 1 for F{z). The details of the separation of variables, including 
existence, uniqueness, and boundary conditions, are given in Ref. [l^ and will be omitted here. We also acknowledge 
a similar, much older, one-dimensional solution given in Ref. Q. The separation of variables solution is 


ri{E, s) = E ^ 


dk A{k)Jd-i{kE^)e sp+iy. 


(61) 


where Jd-i is an ordinary Bessel function of order d — 1, and the amplitudes A{k) are found by taking a Hankel 
transform of the initial ensemble ry(E,0). When the ensemble begins in the microcanonical distribution with energy 
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Eq, we have ri{E,0) = 6{E — Eq), and a closed form expression for A{k) results. The energy distribution r]{E,s), 
subject to r]{E, 0) = S{E — Eq), is then 


v{E,s) 




dfc kJd-iikEQ )Jd-i{kE^)e 


sk'^ 

' SCd+l) 


(62) 


Making use of an identity of Bessel integrals utilized in Eq. (22) of Ref. [l^, we can solve the integral in Eq. (1621) and 
simplify the expression to 


■q{E,s) = 


d + 1 
sE^ 



d — 3 
4 


^d-1 




2(d + l) 




(63) 


where Id-i is a modified Bessel function of order d — 1. Using this energy distribution, we can find the ensemble 
averaged energy as a function of time: 


{E{s)) 


d 

d+lT 


+ \/ Eg S + Eg ■ 


(64) 


Equation (1631) is only valid under the assumption e <C 1. If we begin with an ensemble where e is order unity or 

larger, over sufficiently long time, the slow particles inevitably gain so much energy that the fast particle assumption 

holds and Eq. (I63|) becomes valid asymptotically. We can thus find a universally valid asymptotic energy distribution 

by considering Eq. (EH) or Eq. dMl) in the limit of very large s. Specifically, if fc <C dj^fEg for all ^ 8(d + l)/s, 

1 

which implies that s ^ S^/Eg{d + l)/d, one can approximate Jd-i{kEg) by the lowest order term in its Taylor 
expansion over the non-negligible contributions to the integral in Eq. (I62p . and the solution reduces to 


Va{E,s) 


2ET{d) [ s 






(65) 


where T is the gamma function. One can easily verify that r]a{E,s) is normalized and obeys the Fokker-Planck 
equation. Using the asymptotic energy distribution Eq. (j65p . we find the ensemble averaged energy at a large times 
to be 


{E{s)) 


d 

TTdT’ 


( 66 ) 


The results of this section are summarized as follows. In the quivering limit, correlations in particle energy decay 
over the time scale of a single collision, and as a result, the energy distribution of an ensemble evolves diffusively, 
regardless of the shape and dimensionality of the billiard boundary. Ensembles universally evolve to the asymptotic 
energy distribution given in Eq. (1651) . and ensemble averaged energy asymptotically grows quadratically in time. 
Before discussing the implications and broader context of these results, we comment on the interpretations of the 
coarse grained quantities I and v?. 

If the particular billiard shape is ergodic, then their exists a characteristic ergodic time scale over which ensembles 
uniformly explore the entire billiard boundary. Invoking ergodicity and replacing time averages with phase space 
averages, we deduce that, over time scales greater than the ergodic time scale, 1 will be the billiard’s mean free path, 
and will be the second wall moment M 2 (q) uniformly averaged over the billiard boundary. This implies that, 
over time scales greater than the ergodic time scale, gi and g 2 are time-independent and that a is merely a constant. 
In this case, the expression for gi in Eq. (1521) is equivalent to the wall formula, which was originally used to model 
energy dissipation from collective to microscopic degrees of freedom in nuclear dynamics In non-ergodic billiards, 
or over time scales shorter than the ergodic time scale in ergodic billiards, I and v? will generally be time-dependent 
and can not be interpreted in terms of properties of the billiard shape alone. Nevertheless, they are still well-defined 
properties of the ensemble; I is simply the ensemble’s average free flight distance over the coarse grained time scale, 
and is the average squared wall velocity for the collisions taking place over the coarse grained time scale. 


IV. DISCUSSION 
A. Approximate Quivering 

The quivering limit is most certainly an idealization of time-dependent billiard motion; no real billiard boundary can 
actually move with zero amplitude and period. However, if the idealized system is defined in a physically consistent 
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manner, then we expect that for smaller and smaller a and r, real time-dependent billiards will be better and better 
approximated by quivering billiards. We now clarify how small a and r must actually be for a time-dependent billiard 
to be well-Mproximated by a quivering billiard. 

In Refs. Q and Q, Lieberman, Lichtenberg, and Cohen studied the Fermi-Ulam model numerically and analytically 
using dynamical systems theory. It was shown that the energy evolution of a particle in the Fermi-Ulam model is 
generically diffusive and c an be described by a F okke r-Planck equation for particle speeds such that, using our notation 
from Sec. ins V <C Uc\/L/a. The value Uc\/Lja is associated with the stability of periodic orbits in u-4' space, 
where v and 4* are the particle velocity and wall phase during collisions, respectively. At particle speeds much below 
Uc\/L/a, Refs. Q and Q show that periodic orbits in space are unstable, dynamical correlations are small, and 
trajectories in ?;-4' space are generally chaotic (the language of the day labelled such trajectories stochastic as opposed 
to chaotic). At particle speeds above Uc^/Lja, periodic orbits begin to stabilize, correlations become important, and 
the presence of elliptic islands and invariant spanning curves inhibit energy growth i i- In a one-dimensional 
quivering billiard, correlations vanish, trajectories are stochastic, and particle energy evolves diffusively, so, based on 
Lieberman, Lichtenberg, and Cohen’s work, we see that a quivering billiard is a good description of the Fermi-Ulam 
model when v <C Uc\/Lja. As a becomes smaller and smaller with Uc held fixed, elliptic islands and invariant spanning 
curves move away to regions of larger and larger particle speeds, correlations become smaller and smaller due to the 
more and more erratic wall motion, and quivering becomes a valid approximation for wider and wider ranges of 
particle speeds. As a approaches zero in the idealized limit, the infinitely erratic wall motion destroys correlations, 
elliptic islands and spanning curves occur only at infinite energy, and quivering becomes an exact description for all 
particle speeds. The same reasoning can be applied to higher dimensional time-dependent billiards; as a becomes 
smaller and smaller with Uc held constant, correlations become smaller and smaller and non-diffusive dynamics occur 
at higher and higher energies. We thus claim that when v <C Uc\/hja for all possible particle speeds v that could 
be observed in a simulation or experiment, where Ic is a characteristic free-flight distance, an arbitrary-dimensional 
time-dependent billiard will be approximately a quivering billiard. 

Due to the inevitable increase in particle energy, the speed bound inequality v <C Uc^/lc/a implies that quivering 
will closely approximate a real billiard simulation or experiment only up to some maximum time tmax- The value 
of tmax depends on the particles’ initial energy distribution, but we can estimate its scaling behavior in situations 
where the actual energy distribution is able to evolve the asymptotic distribution given in Eq. (1651) . In such cases, 
the average particle speed at large times can be estimated from the asymptotic ensemble averaged energy given 
by Eq. ((661), and we find v ^ t u^/lc- Substituting this estimate for v into the speed bound inequality yields 
t ^ {Ic/aY''^ (Ic/uc) = We thus have t^ax (Ic/uc) = {IdaY^'^T- As expected, in the quivering 

limit, tmax diverges. 


B. Consistency 


Quivering wall motion corresponds to volume preserving billiard motion with negligible correlations in particles’ 
energy changes. Therefore, if the quivering limit is actually physically meaningful, then the results obtained in Sec. IIIII 
should agree with previous time-dependent billiard literature for the special case of volume preserving billiard motion 
with negligible correlations in energy changes. We now highlight three such examples. 

In Ref. Q, {AE) and are calculated for a single collision in the Eermi-Ulam model, assuming periodic wall 

motion (which corresponds to volume preserving billiard motion on average) and no correlations in the wall velocity 
between collisions. The authors also assume, without explicitly stating, that the wall velocity is an even function 
of time. The expressions obtained in Ref. [J are in fact identical to our expressions for {AEbjt in Eq. (I^Tl) and 
{{AEbY}b, which can be found by adding {AEb}l to Eq. (1^ . under the assumption that all odd moments of the 
wall velocity M 2 n+i vanish. The odd moments vanish in a quivering billiard when we take the quivering limit of wall 
motion defined by an even function of time, so our results agree perfectly with those of Ref. Q . 

Reference [T^ studies the energy evolution of ensembles of independent particles in chaotic adiabatic billiards in 
two and three dimensions. A Fokker-Planck equation to describe the evolution of the energy distribution is proposed, 
and expressions for the corresponding drift and diffusion coefficients are derived. These results are obtained for 
general adiabatic billiard motion, under the assumption that correlations in a particle’s energy changes decay over 
the mixing time scales corresponding to the frozen chaotic billiard shapes. The expressions for gi and g 2 are given in 
terms of a diffusion constant D, and an explicit expression for D is given using the quasilinear approximation - the 
assumption that energy changes between bounces are completely uncorrelated. Under the quasilinear approximation, 
assuming volume preserving billiard motion, the expressions for gi and g 2 in Ref. [l^ are identical to our two and 
three-dimensional expressions for gi and g 2 in Eqs. dsn and dMl), respectively, for ergodic billiards, over time scales 


greater than the ergodic time scale. Our results are thus consistent with those of Ref. [1^. It is remarked in Ref. 


that it is not precisely clear under what conditions the quasilinear approximation will be valid for time-dependent 
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billiards in general, but roughly speaking, the approximation requires the billiard shapes and motion to be “sufficiently 
irregular.” Our results help clarify this issue; the quasilinear approximation is justified when a time-dependent billiard 
is approximately quivering, and the quasilinear approximation is in fact exact, not an approximation, in the quivering 
limit. 

In Ref. @ , it is shown that the velocity distribution for independent particles in a time-dependent irregular container 
is asymptotically universally an exponential. This work assumes an isotropic velocity distribution, volume preserving 
billiard motion, and a three-dimensional billiard. If we assume an isotropic velocity distribution in a quivering billiard, 
we can change variables from energy to velocity in Eq. (I65L and we find the asymptotic velocity distribution /a(v, s) 
in arbitrary dimensions 


/a(v,s) 


I f2{d+l) [rn\‘^ 2(d+i) ,, ,, 

ndT{d) V ^ V^j "" 


(67) 


In agreement with Ref. @ , the isotropic velocity distribution in a quivering billiard is universally an exponential in all 
dimensions. For a three-dimensional chaotic quivering billiard, where s = at and chaotic mixing ensures an isotropic 
velocity distribution, Eq. (1571) is identical to the velocity distribution obtained in Ref. 


C. Fermi acceleration 


Equation (|66l) shows that the ensemble averaged growth is unbounded, increasing quadratically in time. Unbounded 
average energy growth in time-dependent billiards is known as Fermi acceleration. Fermi acceleration was originally 
proposed by Fermi as the mechanism by which cosmic rays gain enormous energies through reflections off of moving 
magnetic fields and since become an active field of research in its own right. The current research generally seeks 
to determine under what conditions time-dependent billiards allow for Fermi acceleration, and to understand how the 
dynamics of sequence of frozen billiard shapes affects the energy growth rate. In Refs. iflll , it was established 
that sufficiently smooth wall motion in the one-dimensional Fermi-Ulam model prohibits Fermi acceleration, and 
that non-smooth wall motion allows for Fermi acceleration that may be much slower than quadratic in time. While 
the one-dimensional billiard is always integrable, higher dimensional billiards allow for integrable, pseudo-integrable, 
chaotic, or mixed dynamics. In Ref. [H, it was conjectured that fully chaotic frozen billiard shapes are a sufficient 
condition for Fermi acceleration in multi-dimensional time-dependent billiards, and the energy growth rate in such 
billiards was thought to be quadratic in time [Sill- It has since been shown that the problem is a bit more subtle; 
certain symmetries in the sequence frozen billiard shapes can prohibit or stunt the quadratic energy growth in chaotic 
billiards 1^. The problem is complicated for non-chaotic multi-dimensional billiards as well. Integrable billiards may 
prohibit [20l [ or allow quad ratic or slower Fermi acceleration, while exponential Fermi acceleration is possible for 
pseudo-integrable billiards 221 and billiards with multiple ergodic components [l^, [H, [H-H^ with possibly mixed or 
pseudo-integrable dynamics. 

Given the complexities observed in the previous literature, our result in Eq. (1661) is surprising; in the quivering limit, 
regardless of the dimensionality or underlying frozen dynamics, time-dependent billiards universally show quadratic 
Fermi acceleration. The apparent contradiction between our work and previous work is due to a difference in the 
limits studied. Both our work and the previous literature, because of the inevitable speed up of particles, analyze 
time-dependent billiards in the adiabatic limit, where the wall speed is much slower than the particle speed. In the 
previous literature, however, the period of billiard oscillations is typically fixed and non-zero (with numerical results 
often presented as a function of the oscillation amplitude), so in the adiabatic limit, the typical time between collisions 
is always much shorter than the billiard’s oscillation period. In our work, the oscillation period approaches zero, so the 
time between collisions is always much larger than the oscillation period, even in the adiabatic limit where particles 
move much faster than walls. 


D. Fixed wall simplifications 


An alternative simplification similar to the quivering billiard has been frequently employed in the literature. The 
so-called static wall approximation (sometimes called the simplified Fermi-Ulam model) was originally introduced in 
3 in order to ease the analytical and numerical study of the Fermi-Ulam model, and through the years has become 
a standard approximation assumed valid for small oscillation amplitudes, often studied entirely in lieu of the exact 
dynamics. See Ref. d, 0, [10 [l^ for example. Using the notation of Sec. m assuming v ^ Uc so that we may 

ignore glancing collisions for the sake of simplicity, the dynamics of the one-dimensional Fermi-Ulam model can be 
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described by the deterministic map, 

Vb = Vb-i - 2u{tb), 

^ ^ , 2L , gih) + 9{tb-i) 

tb — Efc-l H-1-; 

Vb-1 Vb-1 

while the corresponding static wall approximation is given by the deterministic map, 

Vb = Vb-1 - ‘2u{tb), 

2L 

tb — tb-1 H-. 

Vb-1 


(68a) 

(68b) 


(69a) 

(69b) 


In the above maps, Vb-i is the particle’s velocity just before the 6‘^collision, and tb is the time of the 6*^ collision. 
An analogous static wall approximation can be constructed for higher dimensional billiards [l^, [2^, [2^ . Like the 
quivering billiard, the static wall approximation eliminates the implicit equations for the time between collisions by 
holding the billiard boundary fixed. The two models differ because the static wall approximation assumes u(tb) to be 
a well behaved function. It is common practice to consider stochastic versions of the maps (l68ll and (|69l) . where u{tb) 
is replaced hy u{tb + Q for some random variable C d, [H, [l^ m, in, m • The stochastic case simulates the effects 
of external noise on the system and allows one to average over C when determining ensemble averages, which often 
facilitates analytical calculations. 

In Refs, m HI]: Karlis et al. show that the stochastic static wall map and its analogue for the two-dimensional 
Lorentz gas give one half the asymptotic energy growth rate of the stochastic Fermi-Ulam map. This inconsistency 
exists even for small a, so Karlis et al. conclude that (l6^ is not a valid approximation of (1681) . We add that the 
same factor of two discrepancy can be observed between our quivering billiard expression for gi and the corresponding 
expressions obtained from the deterministic static wall maps given in ii [l^ l26j|. In an early study of the Fermi- 
Ulam model, Ref. [3] obtains a drift term that is actually in agreement with the static wall approximation value, but 
a careful reading reveals that the authors make a series of simplifications that inadvertently reduce their Fermi-Ulam 
model to the static wall approximation. Ref. corrects for the energy inconsistency to a high degree of accuracy 
in the stochastic case by introducing the hopping wall approximation. The hopping wall approximation assumes 
wall motion slow enough such that the moving wall’s position at the 6*^ bounce can be approximated by its position 
at the (6 — I)*^ bounce, or by its position at the time of the particle’s collision with the fixed wall just after the 
{b— bounce. This approximation allows g{tb) in Eq. (I68bl) to replaced by either g{tb-i) or g{tb-i + L/vb-i). An 
analogous hopping wall approximation for two dimensions is presented in . Like the static wall approximation, the 
hopping wall approximation eliminates the implicit equations for the time between collisions, which eases numerical 
and analytical study. Based on the hopping wall approximation’s more accurate asymptotic energy growth rate, 
Karlis et al. conclude in Refs. [1^ [2^ that the energy discrepancy between the Fermi-Ulam model and the static wall 
approximation is due to dynamical correlations induced by small changes in the free flight time between collisions 
which are neglected in the static wall approximation. 

Based on the results of this paper, we propose an alternative explanation of the energy discrepancy. The energy 
discrepancy is observed because the static wall approximation is simply unphysical, and it can not accommodate for 
the fact that, due to the relative motion between the particles and walls, collisions with inward moving walls are more 
likely than collisions with outward moving walls. In fact, defining the quivering billiard without the flux factor in the 
biased distribution (so that the biased and unbiased distribution are equal) reproduces the asymptotic energy growth 
rate predicted by the stochastic static wall approximation. Evidently, the last term in Eq. (I68bl) is responsible for 
the bias towards inward moving wall collisions in the exact Fermi-Ulam model, and hopping wall approximation’s 
estimate of this term is responsible for its more accurate energy growth rate. Although the static wall approximation 
is a mathematically well-defined dynamical system, it is an ill-posed physical system for the following reasons. If a 
billiard boundary is truly static such that (j68bl) somehow reduces to (I69b|) . then we must have a —>■ 0. But if a —>■ 0, 
then Uc ^ 0 and the billiard becomes trivially time-independent unless r —>■ 0 as well. However, if both a and r —>■ 0, 
then u{t) can not be a well-behaved function as required by the definition of the static wall map, and, as argued in 
Sec.ini the wall velocity becomes stochastic. This logic seems to be unavoidable; if the walls are to be genuinely fixed, 
then physical consistency demands that the wall motion must be non-existent or stochastic. Based on this reasoning, 
we propose the following conjecture: any physically consistent, non-trivial, fixed wall limit of a time-dependent billiard 
must be physically equivalent to the quivering limit, and the corresponding quivering billiard as defined in this paper 
yields the correct dynamics and energy growth rate (by physically equivalent, we mean equivalent energy and velocity 
statistics). Of particular note, corrections to the free flight time between collisions are not needed to achieve the 
correct energy growth rate. 
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V. EXAMPLES AND NUMERICS 


We now give explicit examples of quivering billiards in one and two dimensions and support the previous sections’ 
analyses with numerical work. Consider first a one dimensional Fermi-Ulam model with one wall oscillating at a 
constant speed. Following the notation of Sec. im the position of the moving wall about its mean position is given by 


git) 


a[-l + 4^'(t)], 0< 4'(t) < i 

a[l_4(vl/(t)_ 1)], i< vl;(i)<l, 


and the corresponding wall velocity is given by 


u{t) = 


Auc, 

-4uc, 


0 < 4'(t) < i 
i < 4'(t) < 1. 


(70) 


(71) 


The numerical analyses of this Fermi-Ulam model are presented in Figs. [3] and SI The histograms in Fig. [3] show of 
the evolution of the energy distribution of 10^ particles of mass m = 1 in a microcanonical ensemble with initial speed 
riQ = 1 at time t = 0, and the curves show the analytical solution for this system in the quivering limit as predicted 
by Eq. (l63l) . For this simulation, we set L = 1.0, a = 10“®, and r = 10“^, which gives Uc = 10“^. We see good 
agreement, with some small deviation apparent beginning at f = 5000. We suspect that the deviation is due to the 
faster particles interacting with the elliptic islands in phase space, which is not accounted for in the quivering billiard. 
By the time t = 15000, a sufficient number of the particles have gained enough energy such that the system is no 
longer approximately quivering. Further energy gain is stunted by elliptic islands, so we see an excess of probability 
(an excess relative to the quivering billiard energy distribution) begin to build up at low energies. Figure S] shows the 
same Fermi-Ulam model, with Uc = 10“^, for successively smaller and smaller values of a and r at time t = 5000. 
As a becomes smaller, we see the actual energy distribution converge to the distribution predicted by the quivering 
billiard. 

The quivering limit of the Fermi-Ulam model given in Eqs. GOl) and (ED is found by following the procedures 
described in Sec. m We first obtain the unbiased distribution, 


P{ub\0) = ^5{ub - 4uc) + ]^5{ub + 4itc), 
and then the biased distribution P(ub\vb-i), 


(72) 


P{Ub\Vb-l) 


1 “ 7^) i^i'^b - 4Uc) + S{ub + 4Mc)] , 

S{ub + 4uc), 


Vb-i > 4 mc 
Vb-l < iUc. 


(73) 


The drift and diffusion terms corresponding to this quivering billiard are found by following the procedures Sec. lIIIDl 
We note that M 2 {q_b) = 16Wc f^e moving wall, and M 2 (qb) = 0 for the stationary wall, so Eq. (1351) yields 
= (1/2) IGrtg. The coarse grained free flight distance is given simply hyl = L, so we find 


giiE) 

g2iE) 


32-\/2m ui 

- -E^ 

L 

64-\/2m ui 3 

- -E^ 

L 


aE^ 


2aE^. 


(74) 



FIG. 3. Energy distribution rj{E, t) of 10® particles following the exact Fermi-Ulam dynamics with small wall oscillation 
amplitude a = 10“® at times t = 100, 1000, 5000, and 15000. The histograms are generated from numerical simulations, and 
the smooth curve is the analytical solution Eq. (1631) for the energy distribution of a particle ensemble in the corresponding 
quivering billiard. 
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FIG. 4. Energy distribution ri{E, t) at t = 5000 of 10® particles following the exact Fermi-Ulam dynamics for successively 
smaller wall oscillation amplitudes a. The histograms are generated from numerical simulations, and the smooth curve is the 
analytical solution Eq. dMli for the energy distribution of a particle ensemble in the corresponding quivering billiard. The case 
for a = 10“® is shown in the t — 5000 plot in Fig. [3] 


8r 8r Sr 8 



(a) t = too (b) t = 1000 (c) t = 5000 (d) t = 15000 


FIG. 5. Energy distribution ri{E, t) of 10® particles at f = 100, 1000, 5000, and 15000 in a quivering billiard corresponding to 
the quivering limit of the Femi-Ulam model used in Eig. [S] The histograms are generated from numerical simulations, and the 
smooth curve is the analytical solution for the energy distribution given by Eq. (pll . 


y 


FIG. 6. The six-circle clover billiard, constructed from sections of six adjacent equi-radii circles. 



The drift and diffusion terms are independent of time, so the rescaled time s is simply s = at. Using the same values 
for L, m, and Uc the we used in the Fermi-Ulam simulation, we find a ~ 4.53 x 10“®. Figure [5] shows the evolving 
energy distribution in the simulated quivering billiard, with the analytical result predicted by Eq. (l63l) superimposed. 
Our analytical solution agrees very well with the numerical simulation. 

For pedagogical purposes, we now construct and simulate a two-dimensional quivering billiard. For the billiard 
shape, we have chosen the six-circle clover introduced in Ref. depicted here in Fig. [51 We set the normal wall 
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velocities along the billiard boundary to be, 


u(q,t) 


Mc|n(q) • x|, 0 < -^{t) < i 
-Uc|n(q) • x|, i < ^'(t) < 1, 


( 75 ) 


where h(q) is the outward unit normal to the wall at q and x is the unit vector in the x-direction. This choice of wall 
velocities gives in the quivering limit, 


P(ub|0,qb) = -6{ub 


Uc|n(qb) • x|) + -S{ub 


■ Mc|n(qfc) ■ x|) 


(76) 


P{ub\vb-i,qb,Ob) 


^b), Vb -1 > Ua\n{qb) • x| 
6{ub + Mc|n(qf,) • x|), Vb-i < ^^^(qf,) • x|. 


(77) 


The six-circle clover constructed from equi-radii circles is fully chaotic 0 , so over time scales greater than the clover’s 
ergodic time scale, is just M 2 {q) averaged uniformly over the billiard boundary. For any q on the boundary, we 
have M 2 (q) = |n(q) • xp, and from Fig. IH we see the outward normals n(q) are distributed uniformly around 

a unit circle, so we have = (1/2) u^. The coarse grained free flight distance I, over time scales greater than the 
ergodic time scale, is just the billiard’s mean free path. For a two dimensional erg odic billiard, the mean free path is 
given by ttA/S, where A is the billiard’s area and S is the billiard’s perimeter If the radius of the circles used to 
construct the six-circle clover is R, then the geometry of Fig. [6] gives A = R^{4:^/3 + tt) and S = AttR. We thus have 
for the drift and diffusion coefficients. 


9iiE) 
52 (P) 


2\/2m itg 

I 

8\/2m 

3l 


1 1 
E 2 = aE 2 



(78) 


where I = i?(-\/3 -I- 7r/4) is the mean free path. 

Figure 0 shows the energy evolution of a microcanonical ensemble of 10® particles in a quivering clover, with the 
distribution Eq. (1631) superimposed. The particles have mass m = 1 and initial energy Eq = 1/2. We constructed the 
clover with circles of radius i? = 1 and set Uc = 6.35 x 10“® to give a ~ 4.53 x 10“®. Again, we see good agreement 
between the distribution predicted by Eq. (1631) and the simulated energy distribution. 



FIG. 7. Energy distribution rj{E, t) of 10® particles at t = 100, 1000, 5000, and 15000 in a two-dimensional chaotic quivering 
billiard. The histograms are generated from numerical simulations, and the smooth curve is the analytical solution for the 
energy distribution given by Eq. (1631) . 


VI. SUMMARY AND CONCLUSIONS 

In this work, we have defined a particular fixed wall limit of time-dependent billiards, the quivering limit, and 
explored the evolution of particles and ensembles in the resulting quivering billiards. We have conjectured that any 
physically consistent, non-trivial, fixed wall limit of a time-dependent billiard must be physically equivalent to the 
quivering limit, and we have shown that the simplifications allowed by a physically consistent fixed wall limit come at a 
price: deterministic billiard dynamics become inherently stochastic. Although quivering is an idealized limit of billiard 
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motion, we have shown that for smaller and smaller oscillation amplitudes and periods, time-dependent billiards 
become better and better approximated by quivering billiards. Billiards that quiver or approximately quiver behave 
universally; particle energy evolves diffusively, particle ensembles achieve a universal asymptotic energy distribution, 
and quadratic Fermi acceleration always occurs, regardless of a billiard’s dimensionality or frozen dynamics. The 
mechanism for this quadratic Fermi acceleration is analogous to a resistive friction-like force, present due to the 
fluctuations induced by the erratic wall motion, as described by the fluctuation-dissipation relation in Eq. ([54)l . 

Through this work, we have gained some insight into issues that have been discussed in the previous literature. 
Namely, we concluded that in the quivering limit, the quasilinear approximation is exact, not an approximation. 
Also, we showed that the often used static wall approximation fails because it is unphysical and can not take into 
account the statistical bias towards inward moving wall collisions. Energy gain in the static wall approximation is a 
purely mixing effect; unbiased fluctuations in particle velocity produce an average increase in particle velocity squared, 
analogous to a Brownian random walk where unbiased fluctuations in position produce an average increase in squared 
distance from the initial position. From this observation, and the fact that the static wall approximation gives one 
half the asymptotic energy growth rate observed in exact systems, we conclude that in the quivering limit, half of the 
average energy gain observed in a time-dependent billiard is due to the mere presence of fluctuations, and half is due 
to the fact that energy gaining fluctuations are more likely than energy losing fluctuations. 

We close by acknowledging that we have not given a rigorous mathematical proof showing that deterministic time- 
dependent billiards become stochastic quivering billiards in the quivering limit. One possible approach toward such a 
proof would be to define some sort of space of time-dependent billiards consisting of systems with different oscillation 
amplitudes and periods, define a metric to give some notion of distance in this space, and prove that particular 
sequences in this space with successively smaller amplitudes and periods are Cauchy sequences. One could then 
determine what properties the space of systems would need to posses in order to assure that these Cauchy sequences 
converge to limits, and then study the limits by studying the sequences that converge to them. Instead of a rigorous 
mathematical approach, we have taken a more intuitive approach and have attempted to justify our work by using 
physical reasoning and by showing consistency with previous results. We hope that the evidence is convincing enough 
to mitigate our mathematical deficiencies. 
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Appendix 



FIG. 8. Incoming and outgoing particle trajectories at the collision location qt in the full and frozen dynamics, assuming a 
collision wall velocity Ub- The full dynamics trajectory is perturbed by an angle 59\ub relative to the frozen dynamics trajectory. 
nt is the outward normal to the boundary at qj. 


B' 



FIG. 9. The geometrical relationship between qt,, q6+i|F, and qf,+i|u6. q6 and qf,+i|F denote the 6*^ and (6+ 1)*^ collision 
locations in the frozen dynamics, respectively, while qt+i \ub denotes the {b + 1)*^ collision location in the full dynamics, ti;, 
and n(,+i are the outward normals to the boundary at q^ and q(,+i|F, respectively. 


Here, we find ||dq6+i|u6||, the magnitude of the perturbation to the frozen dynamics (6+ 1)*^ collision location 
due to the energy gained or lost at the 5*^ collision in the full dynamics. In the frozen dynamics, the collision angle 
db is equal to the angle of reflection. Let 9b + 59\ub be the reflected angle in the full dynamics, assuming a wall 
velocity of Ub at the collision. We denote Vb-i as the incoming particle speed at the 6*^ collision, vt as the velocity 
component tangent to the wall, vp as the reflected particle’s velocity component perpendicular to the wall in the 
frozen dynamics, and vp\ub as the reflected perpendicular velocity component in the full dynamics. The collision 
kinematics give vp\ub = vp — 2ub. The perturbation 59\ub can be found using the geometry in Fig. |S1 Note that 
tan(0b) = ^ and tan(df, + 59\ub) = ■ Expanding ta.n{9b + S9\ub) to first order in 69\ub, we find 


id.w{9b + 59\ub) = 


Vp\Ub 


Vt 
= tan{6b) 


cos^{0b) 


56\ub 


1 


= ^ + 

Vt cos^ (9b) 




(A.l) 
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Noting that vp\ub = vp — 2ub and vp = Vb-i cos 0b, we solve for 59\ub to find 

59\ub = 2cos{9b) ■ (A.2) 

Vb-l 


Figure ini shows the geometry of the 5*^ and (6+1)*^ collisions in both the full and frozen dynamics, where ||dqb+i 
is the length of the line segment C'D'. We assume that 59\ub is small enough such that the wall appears fiat between 
the frozen and full dynamics’ {b + 1)*^ collision locations. The triangle A'B'C in Fig. [^is similar to the triangle 

ABC in Fig. |S1 so we have We note that |A'C"| is the distance between the and (5 + 1)*^ 

collision locations in the frozen dynamics, so we denote |A'(7'| = Lb\F and find 

\B'C'\ = ‘^Lb\F. (A.3) 

Vb-l 


All angles in Fig. [8] can be found in terms of 9b, 06+i|f, and 59\ub. By applying the Law of Sines to the triangle 
B'C'D', we find 


\C'D'\ 


‘2Lb\F 


cos(6>fc) \ub\ 
sin(6»h+i|F) Vb-l' 


We thus have 


ll<5q6+i|Mh|| 


2Lb\F 


cos{9b) \ub\ 

sin(0f,+i|F) Vb-l' 


(A.4) 


(A.5) 
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